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Abstract. - Recent molecular dynamics simulations show that a dilute relativistic gas equili- 
brates to a Jiittner velocity distribution if ensemble velocities are measured simultaneously in 
the observer frame. The analysis of relativistic Brownian motion processes, on the other hand, 
implies that stationary one-particle distributions can differ depending on the underlying time- 
parameterizations. Using molecular dynamics simulations, we demonstrate how this relativistic 
phenomenon can be understood within a deterministic model system. We show that, depending 
on the time-parameterization, one can distinguish different types of 'soft' ergodicity on the level 
of the one-particle distributions. Our analysis further reveals a close connection between time pa- 
rameters and entropy in special relativity. A combination of different time-parameterizations can 
potentially be useful in simulations that combine molecular dynamics algorithms with randomized 
particle creation, annihilation, or decay processes. 
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Introduction. - Understanding the relation between en- 
semble and time averages poses one of the most fundamen- 
tal problems in statistical physics. Ergodicity - the equiv- 
alence of the two averaging procedures - is a commonly 
employed assumption in statistical mechanics [1], albeit 
difficult to prove for realistic systems. During the past 
decades, the ergodicity hypothesis was intensely exam- 
ined for nonrelativistic classical [2-5] and quantum mod- 
els [6-8]. However, much less is known about its meaning 
and validity in relativistic settings [9], when even more 
basic concepts like 'stationarity' may become ambigous as 
time becomes relative [10-12]. A clear conception of the 
interplay between time parameters and thermostatistical 
concepts, like entropy [13-15], is crucial, e.g., if one wishes 
to generalize non-equilibrium fluctuation theorems to a 
relativistic framework [16,17]. Given the rapidly increas- 
ing number of applications in high-energy physics [18, 19] 
and astrophysics [20,21], a firm conceptual foundation is 
desirable not only from a theoretical, but also from a prac- 
tical perspective. 

Ideally, one would like to tackle relativistic many- 
particle problems within a quantum field theory frame- 
work, as this allows for the consistent treatment of parti- 
cle creation, annihilation, or decay processes [22]. In re- 
cent years, substantial progress has been made towards a 



better understanding of both equilibrium [23,24] and non- 
equilibrium processes [25-28] in the context of relativistic 
field theories. However, while without doubt conceptu- 
ally preferable, an exact quantum theoretical treatment is 
in many situations practically unfeasible and, for a con- 
siderable number of applications (e.g., sufficiently dilute 
gases or plasmas), not necessary. With regard to com- 
puter simulations of relativistic systems, suitably adapted 
quasi-classical particle models [29] often provide a more 
efficient basis for quantitive numerical analysis. In partic- 
ular, given the rapid improvement of GPU programming 
tools [30-33] over the past two years, relativistic molecular 
dynamics (MD) simulations as well as other parallelizable 
approaches, e.g., particle based Monte-Carlo (MC) algo- 
rithms [34], can be expected to play an increasingly im- 
portant role in the future. Against this background, the 
present paper addresses selected thermostatistical proper- 
ties of quasi-classical relativistic systems. 

More precisely, we intend to demonstrate that even a 
relatively simple, relativistic model system [35] may pro- 
vide insights into basic conceptual questions, such as: 
How are observer-time and proper-time averages of single- 
particle trajectories related to each other? How are the 
resulting time-averaged distributions linked to stationary 
distributions obtained from simultaneous ensemble mea- 
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surements? Is it possible to establish a connection be- 
tween time paramctcirs and entropy? As wc shall sec, the 
answers also provide further clarification of results that 
were recently obtained in the theory of relativistic Brow- 
nian motions [11,12,36-39]. For example, one can show 
that changing the time-parameterization of a relativistic 
stochastic process (e.g., form coordinate-time to proper- 
time) entails a modification of the corresponding station- 
ary distribution [11]. Below, we will discuss how this phe- 
nomenon can be understood on the basis of a simple de- 
terministic model system. Furthermore, we arc going to il- 
lustrate how thermodynamic state variables become modi- 
fied when replacing coordinate-time through a proper-time 
parameterization. Generally, proper-time paramctcriza- 
tions provide a natural framework for including particle 
decay processes in relativistic simulations, whereas global 
coordinate-time paramctcrizations arc better suited for 
quantifying the many-particle dynamics and, in particu- 
lar, the causal ordering of collision events [29]. Thus, as 
briefly outlined in the latter part of this paper, a combina- 
tion of different time parameterizations may yield useful 
mixed MD/MC-simulation schemes. 

Time- averaged single-particle distributions. Let us 
start by considering the motion of a specific particle in an 
inertial frame Sq. The velocity V := dX/dt of the par- 
ticle can be parameterized in terms of the S-coordinate- 
time t, denoted by V{t), or, alternatively, by the particle's 
proper-time (units such that the speed of light c = 1) 

T(i) = /* diVi - vit'f, (1) 

corresponding to a function V{t) that satisfies V{t) = 
V{T{t)). We may then define the f-averaged velocity prob- 
ability density function (PDF) of the particle by 

ft{v) = jJ^dt'5[v-V{t% (2a) 
and, similarly, the associated r-averaged PDF by 

Mv) = - r dT'S[v-V{T% (2b) 
T Jo 

We would like to understand how the two PDFs ft and 
fr are related to each other as t,T ^ oo. To this end, 
we change the integration variable in Eq. (2b) to the lab 
time, yielding 

M-) = '-i^-r^'M-)='-^^ (3) 

where 7(1;) = (1 — v^)~^/^ is the Lorentz factor. We now 
define stationary distributions 

foc{v) = lim ft{v), foc{v) = lim fr{v). (4a) 

t — ^oc r — ^00 

Analogous to Eq. (3), these equalities are to be understood 



in a distributional sense, i.e., 

J d^v U{v)g{v) = lim J d'^v Mv)g{v), (4b) 

f d^v foo{v)g{v) = lim [ d'^v fr{v)g{v) (4c) 

for any sufficiently well-behaved, physically relevant test 
function g{v). If the dynamics is such that the limits foo 
and foo exist, then Eq. (3) implies that 

/oo(w) = a"^ foo{v)h{v), (5a) 
where the constant 

a = j d^v f^{v)h{v) (5b) 

ensures normalization. Equation (5) states that, asymp- 
totically, the f-avcragcd distribution f^o differs from the 
r-averaged distribution f^o by a factor proportional to the 
relativistic particle energy e = m^{v), where m is the rest 
mass of the particle. Wc shall return to this point when 
discussing the associated entropy functionals further be- 
low. Before doing so, however, let us compare the time- 
averaged PDFs with 'ensemble-averaged' one-particle ve- 
locity PDFs as recently measured in computer experi- 
ments [35,40-42]. 

Ensemble- averaged one-particle distributions. - Molecu- 
lar dynamics simulations by different groups [35, 40-42] 
confirm that the stationary one-particle velocity PDF of 
a d-dimensional dilute relativistic gas in equilibrium is ac- 
curately described by the Jiittner distribution [43,44] 

/j(t;) = Zj-imS(t;)2+'^exp[-/?m7(t^)], \v\<l. (6) 

Here, T = (fce/?)^^ niay be interpreted as a (rest) tem- 
perature, k-g denotes the Boltzmann constant, and Zj = 
Z,]{d,m., 0) the normalization constant. It is worthwhile 
to take a closer look at how exactly the measurements are 
performed in these simulations: 

(i) Velocities are measured in the rest frame Eq of the 
boundary in the case of an enclosed system, or the centcr- 
of-mass frame in the case of periodic boundary conditions. 

(ii) The velocities of all particles are measured t- 
simultancously in Eq, where t is the coordinate-time of Eo- 

This procedure can be interpreted as constructing the 
one-particle velocity PDF via ensemble-averaging ^. In 
the next part, we would like to compare the results of 
this method with those obtained by time- averaging over a 
single-particle trajectory. 

Numerical simulations. - In order to understand how 
ensemble-averaged and time-averaged PDFs are related 
to each other, we consider the fully relativistic (H-l)- 
dinicusioiial two-component gas model studied in [35]. 

-'Simultaneous measurements can be easily made in computer 
sinmlations, but are are very difficult to perform in real experiments 
due to the finiteness of signal speeds in relativity; see discussion 
in [10,42]. 
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In this model, the gas consists of classical, impenetrable 
point-particles (A^i light particles of rest mass mi, and N2 
heavy particles of rest mass m2 > mi). Neighboring parti- 
cles may exchange momentum and energy in elastic binary 
collisions, governed by the relativistic energy-momentum 
conservation laws 

e{mA,PA) + f-[rnB,PB) = e{niA,PA) + <^{niB,PB), 

Pa+Pb = Pa+Pb, (7) 

where A,B {1, 2}, e{m,p) = [in? + p^)^^^, p = mv"/{v), 
and tilde-symbols indicate quantities after the collision. 
Interactions with the boundaries are elastic, i.e., p — > —p 
in the lab frame Eqi defined as the rest frame of the 
boundaries. The dynamics of this system can be ex- 
actly integrated numerically, and the total energy Sq = 
Si=i^(™iiPj) + ^f=i^{''^2,Pj) is conserved in the lab 
frame Sq. We distinguish four types of measurements. 

(a) t-ensemble average: After a period of equilibration, 
we simultaneously measure the velocities of all parti- 
cles in Eq at a given instant of time t. This procedure 
is repeated for energetically equivalent random initial 
conditions mimicking a micro-canonical ensemble at 
energy Sq [35]. 

(b) t-trajectory average: We choose a specific particle of 
either species and measure their velocities at several 
equidistant instants of time t'^^\ . . . , t^"^ 

(c) r-ensemble average: We compute the proper-time 
for each particle during the simulation and measure 
their velocities at a fixed proper-time value ti — . . . ^ 
TN1+N2 = T- Again, this procedure is repeated for 
energetically equivalent random initial conditions. 

(d) T-trajectory average: We choose a specific particle of 
either species and measure their velocities at several 
equidistant instants of proper-time t^^\ ■ ■ ■ , r^"). 

Results. - Figure 1 depicts the equilibrium distributions 
computed from one-dimensional simulations as described 
above. In the case of the ensemble measurements (a) and 
(c) we averaged over 50 different, energetically equivalent 
initial conditions. The single-particle time-averages (b) 
and (c) were determined by measuring velocities at 5 • 10^ 
instants using time intervals At = At = 4-10~'*L/c, where 
L is the system's spatial extension. 

Let us first compare the distribution functions obtained 
by the two t- averaging methods (a) and (b), respectively. 
As evident from the (blue) diamonds and (magenta) cross- 
symbols in Fig. (1), the two different procedures both yield 
a Jiittner distribution with same parameter (3 for either 
species, i.e., 

foo{v)^h{v). (8a) 

Similarly, upon comparing the histograms obtained by 
(c) T-ensemble averaging, see red triangles, and (d) r- 
trajectory averaging, green plus-symbols, we find that 




v/c 



Fig. 1: (Color online) Numerically measured one-particle ve- 
locity PDFs using lab time and proper-time parameterizations, 
respectively. Symbols/methods: (a) t-ensemble average: blue 
o, (b) f-trajectory average: magenta x, (c) r-ensemble average: 
red A, and (d) r-trajectory average: green The results are 
based on simulations with A^i = 5000 light particles of mass 
nil and N2 = 5000 heavy particles with mass 1712 = 2mi. 
The solid curves correspond to Jiittner functions (6) with 
/3 = 0.709 (mic^)~^, but different particle masses, respectively. 
Dashed lines show the corresponding modified Jiittner distri- 
bution, Eq. (8b), with same parameter /3. As the distributions 
are symmetric with respect to the origin, only the positive ve- 
locity axis is shown. 

both methods give the same distribution. But this proper- 
time equilibrium PDF differs from the Jiittner function by 
a factor 1/e, i.e., 

foo{v) = a-'fj{v)/j{v) -: fMj{v). (8b) 

Thus, on the one hand, our simulations confirm the va- 
lidity of Eq. (5) for the one-dimensional two-component 
gas model. One the other hand, Eqs. (8) provide two 
'soft' ergodicity statements on the level of the one-particle 
velocity distributions (we adopt the term 'soft' rather than 
'weak', which is already commonly used in a different con- 
text [8] ) . Evidently, it is necessary to distinguish different 
time-parameters when discussing ergodicity in relativistic 
systems. 

The 'r-stationary' modified Jiittner function (8b) was 
derived earlier in Ref. [45] from a simple collision invari- 
ance criterion. Yet another derivation, based on symmetry 
and entropy arguments, was given in Ref. [14]. However, 
at that time it was not understood that the two distri- 
butions /j and /mj refer to different time parameters, re- 
spectively. In fact, combining the above results with the 
arguments given in [14] reveals an interesting relation be- 
tween time parameters and (relative) entropy in special 
relativity. 

Maximum (relative) entropy principle.- To establish a 
connection between t, t and entropy, let us again gener- 
alize to the case of d-dimensional gas. We shall assume 
that gas consists of N identical particles that are enclosed 
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by a vessel such that their total energy Sq is conserved 
in the rest frame Eo of the vesscil. For eonvenienee, we 
rewrite the velocity distributions /j and /mj in terms of 
the momentum coordinate p = mv^{v) as 



(9a) 
(9b) 



where p° = {m? + p^)^/^ is the relativistic one-particle 
energy, and Zmj = ctZj. Written in this form, it becomes 
obvious that the distributions (9) can be obtained by max- 
imizing the relative entropy functional [14] 



SWs] = - J d''p</.(p)ln 



under the constraints 



N 



= J d'^p 4>{p), 

= j d''pcj>{p)p\ 



(lOa) 

(10b) 
(10c) 



where £s/N is the specific energy mean value at constant 
coordinate-time (s = t) or proper-time (s = r), respec- 
tively (cf. discussion below). The function Ps{p) > in 
Eq. (10a) plays the role of a reference density [46-50] and 
ensures that the argument of the logarithm is dimension- 
less. In 1911, Jiittner [43] obtained the distribution cpj by 
postulating a constant reference density pt{p) = po, cor- 
responding the translation-invariant Lebesgue measure in 
momentum space ^. For comparison, the modified dis- 
tribution (pM.] is obtained by fixing a reference density 
Pt{p) oc 1/p" [14]. It is well-known that 

defines a unique (up to multiplicative factors) Lorentz- 
invariant integration measure in relativistic momentum 
space. Since (f>MJ is related to the Lorentz-invariant time 
parameter t, we may conclude: The symmetry properties 
of the reference density , which appears in the relative 
entropy functional for the stationary distribution, reflect 
the symmetry of the underlying tim.e parameter. 

It is worthwhile to briefly comment on the relation be- 
tween the two different energy "state variables" £t and 
£r in Eqs. (10c). As evident from Eqs. (5) and also il- 
lustrated in Fig. 1, the inverse temperature parameter 
is the same for both the t-stationary distribution (/>j and 
the T-stationary distribution (pMj- If we choose the maxi- 
mum (relative) entropy principle as the starting point for 
deriving and ^mj, the parameter (3 enters as a La- 
grangian multiplier for the energy constraint (10c). Due 
to the (Mastic collision dynamics, the total initial energy £q 



•^Morc precisely, Jiittner [43] considered a Lcbcsguc reference 
measure on relativistic one-particle phase space {{x,p)}; however, 
we can neglect the trivial spatial part of the distribution in our dis- 
cussion. 



is conserved in the lab frame Sq, if one measures tem- 
poral evolution in terms of the coordinate-time t. This 
suggests to identify £t = £o in the case of the constant 
reference density Pt{p) = Po, which yields the f-stationary 
Jiittner function Upon inserting into the energy 
constraint (10c), one then finds that the parameter (3 is 
uniquely determined by 



N 



dp 



\nZ, 



(11a) 



2n 



(d-l)/2 



2m<* ( ) i^(d+i)/2(/?m), (lib) 



where d = 1, 2, 3 is the space dimension and Kn[z) denotes 
the nth modified Bessel function of the second kind [51]. 
Since, for any given initial energy value Eq, the param- 
eter P is fixed by Eq. (11a), it only remains to specify 
the corresponding "proper-time state variable" £r, which 
is to be used in Eq. (10c) when considering a Lorentz- 
invariant reference density Prip) (X l/p° in the relative 
entropy definition (10a). In general, the quantities £r and 
£t differ from each other, because they correspond to av- 
erages at r = const, and t = const., respectively. £t can 
be found by inserting the maximizer of S[(p\Qr], i.e., the 
modified Jiittner distribution 0mj(p;/3), into the energy 
constraint (10b). This leads to^ 



£r 

N 



Zj 



(11c) 



2tt 
Pm 



(d-l)/2 



^(d-i)/2(/3m). (lid) 



Thus, by virtue of Eqs. (11), each energy value £t — £o 
corresponds uniquely to a temperature value T = (/cb/?)"^ 
and a "proper-time energy" 5^. Generally, when using a 
maximum (relative) entropy principle based on reference 
measures that refer to different time parameters, it is im- 
portant to keep in mind that the "left-hand-side" of the 
energy constraint - in our example, the parameter £s in 
Eq. (10c) - must be specified in accordance with the un- 
derlying time parameter s. 



■Hn d = 



2 space dimensions, Eq. (11c) can be rewritten as 
£t = N(m + l/l3) which coincides with the classical non-relativistic 
equipartition theorem (unfortunately, this 'coincidence' holds only 
for d = 2, but not for d = 1,3). 
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Outlook: Applications & extensions.- At this point 
one may well ask whether or not (liffc;rent tinie- 
parameterizations and relativistic MD-type algorithms 
are relevant in practice. Generally, quasi-classical MD 
schemes can be useful for simulating the thcrmalization 
of relativistic gases and plasmas at sufRciently low den- 
sities. In the high-density regime more sophisticated 
methods arc required that take into account quantum ef- 
fects [25-27, 27, 28, 52] in a more detailed manner. For 
space dimensions d > 1, classical particle-particle inter- 
actions cannot be formulated in a fully relativistic way 
anymore [53] - unless one tried to keep track of interac- 
tion fields and retardation effects which is too expensive 
numerically. However, by choosing particle collision crite- 
ria carefully, e.g., by evaluating collisions in the two-body 
center-of-mass frame and using effective cross-sections as 
obtained from quantum theories [22], one can achieve an 
acceptable degree of accuracy in low-to-moderate density 
simulations [42]. 

When applying MD methods to relativistic prob- 
lems, coordinate-time parameterizations provide the more 
suitable framework for the causal ordering of collision 
events [29]. On the other hand, a proper-time parameter- 
ization is the more natural choice if one wishes to model 
decay processes on a classical particle level. This suggests 
to implement hybrid algorithms that combine determinis- 
tic MD schemes for the particle motion with randomized 
particle decay, annihilation, or creation. As a straightfor- 
ward extension of the above algorithm, one could account 
for such events by including chemical reactions that trans- 
form one or more particles into another (set of) particle(s). 
For example, in the simplest case, the proper-time inter- 
vals between successive decay events can be sampled from 
an exponential distribution whose mean life-time is deter- 
mined by quantum field theory. Generally, the creation, 
reaction, and decay rules have to be chosen in agreement 
with the probabilitities and constraints (i.e., conservation 
laws) derived from the corresponding field theories [22]. 

With regard to future investigations, numerical ap- 
proaches appear to be particular promising when com- 
bined with novel GPU programming tools [30-33]. Com- 
pared with conventional CPU-based computations, GPU 
simulations techniques can significantly reduce computa- 
tional costs (up to two orders of magnitude) of paralleliz- 
able MD or MC algorithms. An example of the latter class 
are relativistic MC simulations similar to those performed 
by Peano et al. [34], who generated the relativistic Jiittner 
distribution (6) using randomized collisions. Their algo- 
rithm could be easily generalized to also include quantum 
phenomena such as, e.g., decay processes. In this context, 
it is worth mentioning that the energy distribution of an 
unstable particle at the end of its life-time (assuming the 
particle lives long enough to allow for thcrmalization) will 
be described by (pMj rather than 

Conclusions.- The above discussion shows that, de- 
pending on the application at hand, it may be useful to 
distinguish different types of stationarity in relativistic 



systems. Our numerical results for the one-dimensional 
two-component gas illustrate that the Jiittner distribu- 
tion [43] oc exp(— /3p°) is linked to coordinate-time t, 
whereas the modified distribution (^mj oc exp(— /3p°)/p° 
is linked to proper-time r. They further demonstrate that 
coordinate-time-averaging along a single-particle trajec- 
tory is equivalent to t-simultaneous ensemble averaging 
over many particles. This may be interpreted as a soft 
form of t-ergodicity on the level of the one-particle veloc- 
ity PDF in this model. An analogous statement holds for 
proper-time-parameterizations (soft r-ergodicity) . More- 
over, the deterministic gas model provides a 'micro- 
scopic' illustration of conceptually similar results that 
were recently obtained within the theory of relativistic 
Brownian motions [11]. Similar to deterministic processes, 
relativistic Brownian motions can either be parameterized 
in terms of the coordinate-time t or their proper-time r. 
The resulting stationary momentum distributions (if exist- 
ing) are then connected by a relation equivalent to Eq. (5). 
This result is, perhaps, more difficult to understand (or 
accept) when considering stochastic processes based on 
postulated random driving processes. The deterministic 
model system considered here helps to clarify the micro- 
scopic origin of this relativistic effect. Last but not least, 
our analysis implies a close connection between entropy 
and time parameters, which may be summarized as fol- 
lows: A Lorentz- invariant time parameter corresponds to a 
Lorentz-invariant reference density (probability measure) 
in the entropy functional for the associated stationary dis- 
tribution. 
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